%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all; clc;
%% manages extraction of data on employment bins
dir = "ThirdChapter_revision";
input = ["/Volumes/Transcend/AWFP/",dir,"/data/CalibrationInput/"];
input = join(input,"");
output= ["/Volumes/Transcend/AWFP/",dir,"/data/CalibrationOutput/"];
output= join(output,"");
cd(input)
ende    = 15;
%%
date = {'23062020'}
cd(input)
d       = csvread('delta_avg_estimation.csv');
N       = exp(csvread('count_avg_estimation.csv'));                 
E       = csvread('employment_ts_estimation.csv');
%% eliminate cities for which delta is smaller than 1:
l = linspace(1,size(d,2),size(d,2));
%%
T = size(E,1);
NC= size(E,2);

DD1      = zeros(1,ende);
DD2      = zeros(1,ende);
select8 = zeros(1,ende);
S7 = zeros(ende,NC);

for x = 1:ende

parametersetting;
delta  = d./(1-alpha);                  % transform the firm size pareto tail into the productivity one
varphi = 1 - 1/(gamma*(1-alpha)+1);
Am     = getAm(N,S,alpha,psi,delta);
Edev_t = E - mean(E,1);
rescale= (varphi./Am).*sqrt(D);         % model backed -> use this one

for t = 1:T-1
   Edevresc(t,:) = Edev_t(t,:)./rescale(t,:);   % this is the X
   Edevresc1(t,:)= Edev_t(t+1,:)./rescale(t,:); % this is the Y
end

for i = 1:NC
   rhom(1,i) = (Edevresc(:,i)'*Edevresc1(:,i))/(Edevresc(:,i)'*Edevresc(:,i));
   varrhom(1,i) = 1/size(E,1)*(Edevresc1(:,i)-rhom(1,i)*Edevresc(:,i))'*(Edevresc1(:,i)-rhom(1,i)*Edevresc(:,i));
end


chi  = psi^(1/(1-alpha));
omega= psi^(-1/(1-alpha));

am = (-varrhom + (rhom-1).*(chi+1)+(1-rhom.^2))./((omega-1)*(chi-omega));
cm = (varrhom - (rhom-1).*(omega+1)-(1-rhom.^2))./((chi-1)*(chi-omega));
bm = 1-am-cm;
deltam  = log(am./cm)./log(psi);

select1 = logical(am<1);
select2 = logical(cm<1);
select3 = logical(am>0);
select4 = logical(cm>0);
select5 = logical(bm<1);
select6 = logical(bm>0);
select9 = logical(deltam>1/(1-alpha));
select7 = select1+select2+select3+select4+select5+select6+select9;

cm   = cm(find(select7==7));
am   = am(find(select7==7));
delta   = delta(find(select7==7)); 
deltam  = deltam(find(select7==7));
select8(1,x) = sum(logical(select7==7));
S7(x,:) = select7;

if size(deltam,2)>1
DD2(1,x) = ((1/select8(1,x))*sum((deltam-delta).^2,2))^.5;
aux = corrcoef([deltam',delta']);
DD1(1,x) = aux(2,1);
else
end

end
plot(select8,DD1,'o')
DD2 = [linspace(1,ende,ende)',DD2',DD1',select8']
DD3 = DD2(find(select8>50),:);
%%
clearvars -except date date1 N E d param input output T NC DD3 sel l; clc;
options = optimset('Display','off');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
loss = [DD3,DD3(:,2) + (NC.*ones(size(DD3,1),1)-DD3(:,4))./(NC.*ones(size(DD3,1),1))];
[~,idx] = sort(loss(:,end)); 
sortedloss = loss(idx,:)  
%x = 5       % select here optimal parameter combination
x = 9     
% change here, a bit arbitrary. but you need to motivate the 
% change in solution to be able to use the new (and correct) data
cd(output)
csvwrite('solution.csv',x) 

parametersetting;
delta  = d./(1-alpha);               % transform the firm size pareto tail into the productivity one
varphi = 1 - 1/(gamma*(1-alpha)+1);
Am     = getAm(N,S,alpha,psi,delta);
Edev_t = E - mean(E,1);
rescale= (varphi./Am).*sqrt(D);    % model backed -> use this one

for t = 1:T-1
   Edevresc(t,:) = Edev_t(t,:)./rescale(t,:);   % this is the X
   Edevresc1(t,:)= Edev_t(t+1,:)./rescale(t,:); % this is the Y
end

for i = 1:NC
   rhom(1,i) = (Edevresc(:,i)'*Edevresc1(:,i))/(Edevresc(:,i)'*Edevresc(:,i));
   varrhom(1,i) = 1/size(E,1)*(Edevresc1(:,i)-rhom(1,i)*Edevresc(:,i))'*(Edevresc1(:,i)-rhom(1,i)*Edevresc(:,i));
end

chi  = psi^(1/(1-alpha));
omega= psi^(-1/(1-alpha));

am = (-varrhom + (rhom-1).*(chi+1)+(1-rhom.^2))./((omega-1)*(chi-omega));
cm = (varrhom - (rhom-1).*(omega+1)-(1-rhom.^2))./((chi-1)*(chi-omega));
bm = 1-am-cm;
deltam  = log(am./cm)./log(psi);

select1 = logical(am<1);
select2 = logical(cm<1);
select3 = logical(am>0);
select4 = logical(cm>0);
select5 = logical(bm<1);
select6 = logical(bm>0);
select9 = logical(deltam>1/(1-alpha));
select7 = select1+select2+select3+select4+select5+select6+select9;

cm   = cm(find(select7==7));
am   = am(find(select7==7));
delta   = delta(find(select7==7)); 
deltam  = deltam(find(select7==7));
rhom   = rhom(find(select7==7)); 
varrhom  = varrhom(find(select7==7));
select8 = sum(logical(select7==7));

S = S(find(select7==7));
N = N(find(select7==7));
l = l(find(select7==7));

Enew = [];
for t = 1:size(E,1)
    aux = E(t,:);
    aux = aux(find(select7==7));
    Enew = [Enew;aux];
end
E = Enew;

if size(deltam,2)>1
DD2 =  (1/select8)*sum((deltam-delta).^2,2);
aux = corrcoef([deltam',delta']);
DD1 = aux(2,1);
else
end
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% NEW %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
cd(output)
aux = rescale(:,find(select7==7)).*sqrt(varrhom);
csvwrite('instantvolatilitymodel.csv',aux)
csvwrite('instantvolatilitydata.csv',Edev_t(:,find(select7==7)))
csvwrite('selectcities.csv',select7)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% OLD %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
cd(input)
toexport = [l',zeros(size(am,2),1),deltam'.*(1-alpha),delta'.*(1-alpha),rhom',zeros(size(am,2),1),am',cm',S']; 
cd(output)
csvwrite('toexport.csv',toexport)

aux = deltam'.*(1-alpha);
toexport = toexport(find(aux>1),:);

c1 = mean(toexport,1);
c2 = var(toexport,1).^.5;
c3 = min(toexport,[],1);
c4 = max(toexport,[],1);

table = [c1(:,2:end);c2(:,2:end);c3(:,2:end);c4(:,2:end)];
csvwrite('tablemoments.csv',table)

table1 = [prctile(S,10),prctile(S,50),prctile(S,90);prctile(am,10),prctile(am,50),prctile(am,90);prctile(bm,10),prctile(bm,50),prctile(bm,90);prctile(cm,10),prctile(cm,50),prctile(cm,90)];
tline = linspace(1,T-1,T-1);
%% 
clear Am Enew Edev_t Edevresc Edevresc1 Emean tline rescale select1 select2 
clear select3 select4 select5 select6 select7 select8 select9 chi omega  
clear D DD1 DD2 toexport bm d

A1 = getAm(N,S,alpha,psi,deltam)./N;
A3 = getAm(N,S,alpha,psi,log(mean(am)./mean(cm))./log(psi))./N;    % variable S
A4 = getAm(N,mean(S),alpha,psi,log(am./mean(cm))./log(psi))./N;    % variable am
A5 = getAm(N,mean(S),alpha,psi,log(mean(am)./cm)./log(psi))./N;    % variable cm

cd(output)
toexport1 = [A1',A3',A4',A5'];
csvwrite('toexport1.csv',toexport1)

